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la  previous  investigations,  a  sequential  estimation  procedure  for  estimating  the 
State  of  a  lunar  orbiting  space  vehicle  acted  upon  by  unmodeled  terms  in  the  lunar 
potential  and  unmodeled  mechanical  forces  such  as  those  due  to  venting  has  been 
developed.  The  unmodeled  accelerations  are  approximated  by  a  first-order  Gauss- 
Markov  process.  The  estimation  algorithm  gives  an  estimate  of  the  position  and 
velocity  of  the  space  vehicle  as  well  as  the  components  of  the  unmodeled  acceleratioi 
at  each  observation  epoch.  The  results  obtained  by  processing  real-time  tracking 
data  from  the  Apollo  10,  11,  and  15  missions  have  indicated  that  the  algorithm 
provides  more  precise  estimates  of  the  vdi  icle  state  than  conventional  orbit 
determination  procedures  and,  hence,  provides  accurate  input  for  navigation  or 
guidance  purposes.  The  question  of  the  usefulness  of  the  unmodeled  accelerations 
for  scientific  experiments  has  not  been  established.  This  investigation  considers 
the  accuracy  with  which  the  algorithm  can  estimate  the  acceleration  due  to  modeled 
lunar  surface  mascons  by  numerical  simulation  of  the  estimation  process.  It  is 
shown  that  an  accurate  estimate  of  the  history  of  the  unmodeled  acceleration  can  be 
obtained.  The  investigation  also  considers  the  effects  of  the  magnitude  and  location 
of  the  mascons,  as  well  as  the  effect  of  the  observation  accuracy. 
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Abstract 


In  previous  investigations,  a  sequential  estimation  procedure  for  estimating  the 
state  of  a  lunar  orbiting  space  vehicle  acted  upon  by  unmodeled  terms  in  the 
lunar  potential  and  unmodeled  mechanical  forces  such  as  those  due  to  venting  has 
been  developed.  The  unmodeled  accelerations  are  approximated  by  a  first-order. 
Gauss-Markov  process.  The  estimation  algorithm  gives  an  estimate  of  the  position 
and  velocity  of  the  space  vehicle  as  well  as  the  components  of  the  unmodeled  ac¬ 
celeration  at  each  observation  epoch.  The  results  obtained  by  processing  real¬ 
time  tracking  data  from  the  Apollo  10,  11,  and  15  missions  have  indicated  that 
the  algorithm  provides  more  precise  estimates  of  the  vehicle  state  than  con¬ 
ventional  orbit  determination  procedures  and,  hence,  provides  accurate  input  for 
navigation  or  guidance  purposes.  The  question  of  the  usefulness  of  the  unmodeled 
accelerations  for  scientific  experiments  has  not  been  established. 

This  investigation  considers  the  accuracy  with  which  the  algorithm  can  estimate 
the  acceleration  due  to  modeled  lunar  surface  mascons  by  numerical  simulation 
of  the  estimation  process.  '  It  is  shows  that  an  accurate  estimate  of  the  history 
of  the  unmodeled  acceleration  can  be  obtained.  The  investigation  also  considers 
the  effects  of  the  magnitude  and  location  of  the  mascons,  as  well  as  the  efftct 
of  the  observation  accuracy. 
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Introduction 


The  problem  of  estimating  the  state  of  a  nonlinear  dynamical  system  in¬ 
fluenced  by  random  forces  using  discrete  observations  which  are  subject  to  random 
error  has  received  considerable  interest  during  the  past  decade.  Impetus  for 
this  interest  has  been  derived  from  requirements  for  precise  determination  of  the 
orbits  and  trajectories  of  space  vehicles.  The  precise  determination  of  the  mo¬ 
tion  of  a  space  vehicle  is  necessary,  first,  so  that  the  navigation  functions 
required  to  establish  a  specified  trajectory  can  be  performed  and,  second,  so 
that  measurements  of  the  satellites'  motion  can  be  used  to  infer  properties  of 
the  physical  medium  in  which  the  vehicle  is  moving.  The  requirements  of  deter¬ 
mining  the  motion  of  a  lunar  orbiting  space  vehicle  have  placed  extensive  demands 
on  the  conventional  least  squares  orbit  determination  method.  In  particular,  the 
problem  of  determining,  simultaneously,  precise  estimates  of  the  vehicle  posi¬ 
tion  and  velocity  while  estimating  the  values  for  the  coefficients  of  the  lunar 
gravity  field  may  require  accuracies  which  exceed  the  capabilities  of  classical 
orbit  determination  methods. 

t  , 

The  solution  to  the  classical  orbit  determination  problem  involves  linear¬ 
izing  the  nonlinear  equations,  which  define  the  problem,  about  a  specified 
reference  trajectory  and  then  applying  linear  estimation  techniques.  Cnee  the 
problem  has  been  reduced  to  a  linear  estimation  problem,  the  estimation  algorithm 
is  obtained  by  applying  the  weighted  least  squares  method  or  by  use  of  statistical 
based  estimation  methods,  such  as  minimum  variance  or  maximum  likelihood  methods. 
Using  the  results  from  either  of  these  approaches,  the  estimate  may  be  obtained 
by  using  a  batch  data  processing  algorithm,  in  which  all  of  the  observations  are 
processed  simultaneously  to  obtain  an  estimate  of  the  state  at  some  reference 
epoch,  or  the  state  estimate  may  be  obtained  by  processing  each  data  point  sequen¬ 
tially  as  it  is  obtained. 


4 


It  is  well  known  that  errors  of  four  basic  types  influence  the  accuracy  of 
the  estimate  obtained  by  each  of' these  procedures;  i.e.,  1)  errors  due  to  the 
linearization  assumptions,  2)  errors  introduced  in  the  computational  procedure, 

3)  errors  which  occur  in  the  observation  process,  and  4)  errors  due  to  inaccuracies 
in  the  mathematical  model  used  to  define  the  problem.  The  effect  of  these  errors 
on  batch-type  estimation  algorithms  is  discussed  in  Ref.  (1),  (2),  and  (3).  If 
the  extended  form  of  the  sequential  estimation  algorithm^ 4 ^  is  used,  where  the 
estimate  of  the  state  at  a  time,  t,.  ,  is  used  to  define  the  reference  trajectory 
for  propagating  the  estimate  from  the  observation  epoch,  t..  ,  to  the  next  ob¬ 
servation  epoch,  t j ,  the  effects  of  the  nonlinearities  are  minimized.  Further¬ 
more,  in  most  orbit  determination  problems,  sufficient  accuracy  can  be  obtained 
during  the  computation  process  to  eliminate  the  integration  errors.  Hence,  using 
the  extended  sequential  estimation  algorithm,  the  effects  of  the  last  two  error 
sources  will  be  the  primary  factors  which  limit  the  orbit  determination  accuracy 
for  near-earth  or  near-lunar  satellites.  While  the  effects  of  the  observation 

errors  are  an  important  factor,  the  increase  in  the  accuracy  of  the  doppler 

•  , 

tracking  data  coupled  with  the  high  precision  laser  ranging  data  suggests  that 
errors  in  the  data  are  less  important  in  limiting  the  accuracy  obtained  with  con¬ 
ventional  orbit  determination  procedures  than  are  the  errors  due  to  model  inac¬ 
curacies.  As  a  consequence  of  this  fact,  primary  attention  in  this  investigation 
will  be  directed  towards  consideration  of  methods  for  reducing  the  errors  due  to 
inaccuracies  in  the  dynamic  model. 

It  should  be  noted  that  there  is  some  trade-off  between  the  errors  which  oc¬ 
cur  in  the  computational  process  and  the  errors  due  to  an  incomplete  or  inaccurate 
mathematical  model.  With  a  more  complete  mathematical  model,  it  is  more  difficult 
to  obtain  an  accurate  numerical  solution  to  the  relations  which  define  the  estim¬ 


ation  procedure.  On  the  other  hand,  if  the  mathematical  model  is  simplified  to 
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alleviate  the  computational  problems,  important  physical  effects  may  be  neglected. 

The  subsequent  estimation  procedure  may  "diverge"  due  to  model  error^1^5).  This 

difficulty  may  be  encountered,  also,  if  the  actual  dynamical  process  is  not 

understood  well  enough  to  allow  a  precise  definition  of  the  mathematical  model. 

The  effects  of  inaccuracies  in  the  mathematical  model  on  the  accuracy  and  the 

stability  of  the  linear  sequential  estimation  procedures  are  discussed  in  Ref.  (4), 

(5),  and  (6).  A  number  of  methods  have  been  proposed  to  alleviate  the  problem  of 

estimation  divergence.  These  methods  include  additions  to  the  state  error  co- 

variance  matrix  to  account  for  noise  in  the  differential  equations  which  govern 

the  motion^ specification  of  a  minimum  bound  on  the  estimation  error  covariance 

matrix^8),  and  the  utilization  of  a  finite^ **)  or  decaying^ 8 ^  data  set.  While 

these  methods  lead  to  a  more  stable  estimation  algorithm,  they  suffer  from  the 

disadvantage  that  the  accuracy  of  the  estimate  is  determined  by  certain  empirical 

parameters  which  must  be  specified  "apriori"  without  knowledge  of  the  disturbing 

process.  In  addition,  the  methods  do  not  yield  information  directly  related  to 

the  unmodeled  accelerations.  Such  information  is  of  considerable  value  in  post- 

»  . 

flight  data  analysis  for  improving  the  knowledge  of  the  mathematical  model. 

In  the  following  discussion,  a  sequential  estimation  method  which  compen¬ 
sates  for  unmodeled  effects  in  the  differential  equations  of  the  dynamical  process, 
is  described.  The  method  has  two  advantages:  1)  it  can  be  used  to  obtain  an 
improved  estimate  of  the  vehicle  state  during  real  time  estimation  problems  and 
2)  the  method  will  yield  information  which  can  be  used  in  post- flight  data 
analysis  to  improve  the  mathematical  model.  In  the  proposed  method,  the  "un¬ 
modeled"  accelerations  are  assumed  to  consist  of  the  superposition  of  a  time  cor¬ 
related  component  and  a  purely  random  component  and  are  approximated  by  a  first- 
order  Gauss-Markov  process.  This  model  has  been  used  previously  to  compensate  for 


__ • _ irm rMatl'mi^liin  , 
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errors  in  inertial  navigation  devices and  to  compensate  for  t're  effects  of 
random  fluctuations  in  atmospheric  drag^ 1 1  ^ . 

In  Refs.  (10)  and  (11),  a  stationary  Gauss-Karkov  process  is  assumed  and 
the  correlation  time,  or  equivalently,  the  correlation  coefficient  is  assumed  to 
be  constant  and  known,  apriori.  In  Refs.  (12)  and  (13),  this  assumption  is 
relaxed  and  the  correlation  coefficients  are  regarded  as  unknown  parameters  and 
the  best  estimate  of  their  values  are  determined  during  the  estimation  process 
In  Refs.  (12)  and  (13),  this  approach  is  used  to  develop  a  sequential  estimation 
procedure  for  estimating  the  state  of  a  lunar  orbiting  space  vehicle  acted  upon 


i 


by  unmodeled  forces  due  to  an  incomplete  lunar  potential  and  unmodeled  mechanical 
forces  such  as  those  due  to  venting,  water  dumps,  or  translational  forces  due  to 
unbalanced  attitude  control  reactions.  The  estimation  algorithm  gives  an  estimate 


of  the  position  and  velocity  of  the  space  vehicle  as  well  as  the  componen  :s  of 

the  unmodeled  acceleration  at  each  observation  epoch.  The  algorithm,  referred 

to  as  the  Dynamic  Model  Compensation  (DKC)  method,  has  been  applied  in  further 

simulated  studies  in  Refs.  (14)  and  (15).  The  results  of  these  studies  along 

•  . 

with  the  results  from  processing  real  time  tracking  data  from  the  Apollo  10,  11, 
and  15  missions^ 13^  have  indicated  that  the  algorithm  provides  more  precise 


estimates  of  the  vehicle  state  than  the  conventional  least  squares  batch  or 
sequential  orbit  determination  procedures  and,  hence,  provides  a  more  accurate 


input  for  navigation  or  guidance  purposes.  The  results  obtained  in  these 
studies  indicate  that  the  estimation  procedure  is  stable  in  the  sense  that  the 
error  norm  does  not  grow  in  an  unbounded  manner  and  that  the  norm  of  the  covariance 
matrix  provides  a  conservative  bound  on  the  error  in  the  estimate.  Finally,  the 
estimates  of  the  unmodeled  accelerations  obtained  by  processing  range-rate  tracking 
data  from  the  lunar  orbit  phase  of  the  Apollo  10  and  11  missions  indicated  a  nigh 
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correlation  with  the  location  of  lunar  surface  mass  concentrations ,  as  reported 
in  Refs.  (12^  and  (13).  In  these  studies  Lhe  observation  residual  was  bounded 
by  the  formal  standard  deviation  assigned  to  the  range-rate  measurements  and, 
consequently,  it  is  conjectured  that  the  estimates  of  the  unmodeled  accelera¬ 
tions  are  reasonable  representations  of  the  true  unmodeled  accelerations  acting 
on  the  vehicle.  However,  since  actual  tracking  data  is  used  for  the  studies 
reported  in  Refs.  (12)  and  (13),  the  true  acceleration  history  is  not  available 
for  comparison. 

The  primary  objective  of  the  investigation  described  in  the  following 
sections  is  to  determine  the  accuracy  with  which  the  unmodeled  accelerations  can 
be  estimated  using  the  Dynamic  Model  Compensation  (DMC)  method  proposed  in 
Refs.  (12)  and  (13).  The  question  is  answered  by  simulating  numerically  the 
orbit  determination  procedure  for  a  space  vehicle  in  an  Apollo-type  lunar  orbit. 

In  the  simulated  study,  consideration  is  given  to  the  effects  of  the  accuracy  and 
density  of  the  observations,  as  well  as  the  magnitxide  of  the  unmodeled  accelerations, 
on  the  accuracy  of  the  estimate. 

t  . 
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The  Estimation  Algorithm 

The  equations  which  describe  the  notion  of  a  lunar  satellite  can  be  ex¬ 
pressed  by  the  following  system  of  first-order  differential  equations: 

r  =  v  ,  v  =  a  +  a  +m  (1) 

c  p 

where  r  is  a  three-vector  of  selenocentric  position  components,  v  is  a  three- 
vector  of  selenocentric  velocity  components,  ac  is  the  acceleration  due  to  the 
central  body,  a^  is  the  modeled  acceleration  due  to  other  sources,  gravitational 
or  otherwise,  and,  the  three- vector ,  m  ,  represents  the  effects  of  all  accelera¬ 
tions  not  accounted  for  in  the  mathematical  model  used  to  describe  the  motion  of 
the  satellite.  In  the  subsequent  consideration,  m  will  be  referred  tr  as  the 
"unmodeled  acceleration". 

In  the  Dynamic  Model  Compensation  procedure  described  in  the  following 
discussion,  the  unmodeled  acceleration  m  is  approximated  as  an  adaptive, 
first-order,  Gauss-Markov  process,  e(t)  ,  which  satisfies  the  following  vector 
differential  equation  , 


e(t)  =  Be(t)  +  u(t) 


where  e(t)  is  a  three-vector  whose  components  approximate  the  values  of  the 
unmodeled  accelerations  at  the  time,  t  ,  and  u(t)  is  a  three-vector  of 
Gaussian  noise  whose  components  are  assumed  to  satisfy  the  apriori  statistics: 


E[u(t)]  =  0  ,  E[u(t)u  (t)]  -  q(t )<5(t  -  t) 


where  <5(t  -  t)  is  the  Dirac  "Delta"  function  and  where  q(t)  is  assumed  to  oe 
a  diagonal  matrix.  This  assumption  implies  that  the  components  of  u  are  not 
correlated.  The  coefficient  matrix,  B  ,  is  defined  by  the  components  = 

(i,  j  =  1,  2,  3)  where  the  8^  are  assumed  to  be  unknown  parameters  whose 


... 
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values  are  to  be  determined  during  the  estimation  process  by  the  inclusion  of  the 
equation  6  =  0*  where  8  -  C ^  e  substituted  for  m  in  Eq.  (1) 

and  if  the  state  vector  X  is  defined  as 

„T  r  T  T  T  eT,  / 1.  \ 

X  =  [r  :  v  :  e  :  8  J  1*0 

then  Eqs.  (1),  (2),  and  the  relations  8=0  can  be  combined  to  obtain  the  fol¬ 
lowing  differential  equations  for  the 'state  of  the  dynamical  system: 

X  =  F(X,  u,  t)  ,  X(tQ)  =  XQ  (5) 

T  T  T  T 

where  F  =  [v  :  (a  +  e)  :  (Be  +  u)  :  0]  and  where  a  =  a  t  a  ,  i.e.,  the 

m  m  c  p 

components  of  the  modeled  acceleration.  In  the  usual  orbit  determination  problem. 


the  initial  conditions ,  X  ,  are  unknown . 

o 

For  t  >  t^  ,  where  t^  is  some  reference  epoch,  the  solution  to  Eq.  (5) 
can  be  expressed  in  integral  form  as  follows: 


=  r.  +  v.it  t 

1  x  Jt. 


a(x)[t  -  x]dx 


=  v.At  +  [ 

1  Jt. 


a(x)dx 


e(t)  =  E(t)ci  +  ,  8(t)  =  8i 


where  At  =  t  -  t.  and  a(t)  =  a  (t)  +  c(t)  .  The  matrices  E(t)  and  l,  are 
l  m  i 

defined  as 


a  0  0 

x 

E(t)  =  0  a  0 

y 

0  0  a 
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t! 

l 


a  *4-a2  u  * 

y  y  y  • 


a 

z 


Sl-a*  u  ] 
z  z 


(8) 


where  a  =  exp  [-8  (t  -  t.)3  with  similar  definitions  for  a  and  a  . 

X  *  X  i  y  2 

Further  details  on  the  development  of  the  third  of  Eqs.  (6)  are  given  in  Ref.  (16). 

Using  the  definition  in  Eq.  (4),  the  solution  to  Eqs.  (5)  can  be  expressed 
as 


X(t)  =  6(X^,  t^,  t)  +  n.  ,  t  >  t^ 


(9) 


T  T  T  T 

where  n.  =  Cn  :  n  :  n  :  03.  The  components  of  the  state  noise  matrix,  n.  , 
i  r  v  c  l 

are  due  to  the  purely  random  components  of  the  unmodeled  accelerations  and  can  be 
defined  as  follows: 

ft 


=  Mt  -  tMt  ,  nv  =  |t 


dr  ,  n  =  l. 
e  l 


(10) 


In  view  of  Eqs.  (3),  the  random  process  n^.  will  satisfy  the  conditions 


E[n.3  =  0  ,  ECmnp  =  Ql  6 

3  t  ,  * 


(11) 


where  5^  is  the  Kronecker  delta.  Eq.  (9)  is  used  to  propagate  the  state  vector 
from  an  observation  point  t^  to  an  observation  point  t_.  . 

The  relationship  between  the  p-dinensional  observation  vector  Y.  ,  the 
p-vector  of  observation  noise,  v.  ,  and  the  state  at  the  time,  t.  ,  is 

r  i  1 


Y.  =  G(X. ,  t.)  +  v.  . 

i  l  i  l 


(12) 


In  the  following  discussion,  it  is  assumed  that  the  observation  noise,  , 
satisfies  the  following  conditions: 


E[v. ]  =  0  ,  ECv.vT]  =  R.6..  ,  E[v.xT]  =  O' 

i  i  5  i  i]  i] 


(13) 
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where  6. .  is  the  Kronecker  delta. 

*D 

The  problem  considered  then  is  posed  as  follows :  Given  the  relation  for 
propagating  the  state,  Eq.  (9),  the  observation-state  relation,  Eq.  (12),  the 
sequence  of  observations  ,  i  =  1,  ...,  k,  and  the  statistics  on  the  state  noise, 
Eq.  (3),  and  the  observation  noise,  Eq.  (13),  find  the  best  estimate,  in  the  min- 

A 

’mum  variance  sense,  of  the  state,  ,  at  the  time  t^. 

Under  the  conditions  given  in  the  previous  problem,  the  estimate  at  the 
time,  t^  ,  caji  be  obtained  using  the  following  algorithm^4 ) : 


V,(WVi'  V 

\  =  4<tV  'k-i'pk-;  4  (t;<’  tk-i)  +  9k-i 

•s, 

A 


k 

vICHk  K  <  *  \r‘ 


P,  = 


*k  +  •Sc  CYk  -  G(*k’  V] 
CI  -  “k  V  pk 


(14) 


"here  =  DG/SXj.]*,  ^  =  ECX^  |  Y^  ....  y]  and  ^  =  ECX,.  |  Y^  ....  Y^]. 
The  covariance  matrices  ?k  and  are  associated  with  the  state  estimates 

A 

and  X^  ,  respectively.  The  state  transition  matrix,  $(tk,  ^)  satisfies 

the  following  differential  equation 


$(t,  tk)  =  A(t)  4(t,  tk)  «(tk,  tk)  =  I 


(15) 


where  A(t)  =  [3F/9X]*  .  The  symbol  [  ]*  indicates  that  the  quantity  in  the 

A 

bracket  is  evaluated  on  the  solution,  X(t)  =  0(^K  ^  q*  ^  <  The 

A 

solution  9(Xk  2*  ^k  l*  ^  Stained  Sq.  (9)  using  the  condition  that 

i  -  E[£  |  Y  ,  ... ,  Y  ]  =  0  .  The  algorithm  given  by  Eqs.  (14)  and  (15)  is  the 
extended  form  of  the  usual  linear  sequential  estimator  or  the  Kalman-3ucy  filter 


g  ' 


12 


as  discussed  in  Ref.  (4).  The  covariance  matrix  is  associated  with  the  best 


estimate  of  based  on  k-observations  while  is  the  covariance  matrix  as¬ 

sociated  with  the  best  estimate  of  X^  based  on  (k-l)-observations. 

* 

It  should  be  recalled  that  the  estimate  X^  ,  includes  an  estimate  at  the 


time,  t^  ,  of  the  components  of  the  position,  the  velocity,  the  unmodeled  ac¬ 


celeration,  e  ,  and  the  correlation  coefficients,  8^,  8^.,  and  8 z«  The  algorithm 


requires  apriori  or  "initial"  estimates  of  each  of  these  quantities  as  well  as 


the  apriori  covariance  matrices,  P  ,  Q. ,  and  R. ,  associated  with  the  initial 

0  1  1 


state  error  and  with  the  state  noise  and  the  observation  noise,  respectively. 

The  development  of  the  algorithm  and  the  computational  procedure  required  to  im¬ 
plement  the  algorithm  are  discussed  in  greater  detail  in  Refs.  (12)  and  (13). 

In  addition  to  processing  the  tracking  data  from  the  Apollo  10  and  11  mis¬ 
sions,  the  investigation  described  in  Refs.  (12)  and  (13)  evaluate  the  char¬ 
acteristics  of  the  estimation  algorithm  using  simulated  data.  It  is  concluded 
from  the  simulated  studies  that  in  the  presence  of  the  unmodeled  accelerations , 

the  algorithm  is  stable  and  the  covariance  matrix  provides  a  reasonable,  but  con- 

»  , 

servative  estimate  of  the  actual  error.  Furthermore ,  the  estimate  of  the  posi¬ 
tion  and  velocity  is  more  accurate  than  the  estimate  obtained  using  the  usual 
sequential  estimator  or  Kalman  Filter,  modified  as  discussed  in  Ref.  (4),  to 
account  for  the  effects  of  state  process  noise.  In  the  simulated  study  which  fol¬ 
lows,  the  accuracy  with  which  the  unmodeled  acceleration  components  can  be  es¬ 
timated  by  the  DKC  procedure  is  investigated. 


. . .  . . . . » - - - - 
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Generation  of  the  Simulated  Observations 

The  simulated  observations  are  generated  by  the  integration  of  the  system 
of  equations 

•T  =  VT  •  *T  *  am  +  "  (16) 

where  r^Ctg),  and  m(t)  are  specified  and  where  in  the  simulated  study 

represents  a  three-vector  of  the  true  position  components,  is  a  three- 

vector  of  the  true  velocity  components,  a  represents  the  terms  modeled  in  the 

estimation  algorithm,  and  m  represents  the  terms  not  directly  included  in  the 

dynamic  model  used  in  the  estimation  algorithm.  The  term  m  will  define  the  true 

model  error  and  is  approximated  by  a  first-order  Gauss-Harkov  process  in  the  DMC 

estimation  algorithm.  The  model  error  is  simulated  in  the  following  study  by 

assuming  that  m  is  the  acceleration  due  to  several  point  masses  at  various 

locations  on  the  lunar  surfa  ■>.  In  the  simulation,  the  term  a  includes  the 

m 

central  body  acceleration  of  the  moon  and  the  acceleration  due  to  a  point-mass 
earth.  The  effects  of  the  sun  and  the  other,  planets  were  neglected. 

The  numerical  integration  of  Eq.  (16)  yields  a  solution  for  the  "true" 
motion  of  the  satellite  which  is  used  in  generating  simulated  observations.  The 
computer  program  ascertains  which  tracking  stations  can  observe  the  satellite 
and  a  range-rate  observation,  p,  is  generated  using  the  position  and  velocity  of 
the  satellite  relative  to  a  topocentric  earth-fixed  tracking  station.  The 
"true"  value  of  the  range-rate  observation  is  corrupted  to  simulate  the  actual 
observation  by  adding  the  quantity  u^X  ,  where  a ^  is  the  standard  deviation 
of  the  observation  noise  and  X  is  a  Gaussian  distributed  random  variable  with 
zero  mean  and  unit  variance. 


The  numerical  integration  used  in  both  the  simulations  and  the  estimation 
procedure  is  a  variable  step  Runge-Kutta  method  formulate^  in  Ref.  (18).  As 
implemented,  integrators  of  two  different  orders  can  be  used.  During  periods  in 

which  observations  are  generated,  a  fixed  step  size  is  used  to  simulate  a  fixed 
observation  interval.  This  interval  is  normally  sufficiently  small  to  allow  use 
of  a  lower  order  integrator,  namely,  the  Runge-Kutta-Fehlberg  4(5).  If  the  ob¬ 
servation  interval  results  in  a  step  size  which  is  too  large  and  the  4(5)  is 
inadequate  for  maintaining  accuracy,  a  higher  order  7(8)  method  car;  be  used. 

The  estimation  program  requires  the  numerical  integration  of  Eqs.  (5)  and 
(15),  whereas  the  observation  generation  requires  the  integration  of  Eq.  (16),  only. 
The  state  vector  for  Eq.  (5)  consists  of  12  elements,  while  Eq.  (16)  involves  six 
elements.  The  state  transition  matrix,  defined  in  Eq.  (15),  consists  of  144 
elements;  however,  the  order  of  the  system  can  be  reduced  by  noting  that  a  number 
of  elements  can  be  integrated  analytically. 

To  simulate  the  real-world,  the  nominal  or  reference  orbit  initial  con¬ 
ditions  were  assumed  to  be  different  from  those  used  in  the  generation  of  the 

^  i  •  •  •  •  • 

observations.  The  differences  between  the  two  sets  of  initial  conditions  for 
all  orbits  considered,  were  as  follows: 

Ax  =  -305.8  meters  ,  A k  =  0.0542  meters/sec 

Ay  =  -304.7  meters  ,  Ay  =  -0.3682  meters/sec 

Az  =  -  71.9  meters  ,  Az  =  0.3135  meters/sec 

Two  orbit  inclinations  were  considered,  i.e.,  I  =  180°  and  I  =  150°.  For  both 
cases,  the  satellite  starts  at  approximately  90°  East  selenographic  longitude, 
with  altitude  of  100  km,  and  with  an  eccentricity  of  zero.  For  the  case  I  =  180°, 
the  satellite  will  be  in  a  lunar  equatorial  orbit.  The  lunar  groundtrack  of  the 


inclined  orbit  is  shown  in  Fig.  1.  Fig.  1  also  shows  the  locations  of  four  mascons 

whose  effect  was  included  in  the' generation  of  the  simulated  observations,  but  was 

not  included  in  the  dynamic  model,  a  ,  for  the  estimation  algorithm.  The  figure 

m 

on  v.'hich  the  simulated  groundtrack  is  drawn  (Fig.  1)  was  provided  by  W.  L. 

Sjogren  of  the  Jet  Propulsion  Laboratory18. 

The  simulated  mascons  shown  in  Fig.  1  are  given  in  Ref.  (18).  (See  also 
Table  I).  It  should  be  noted  that  the  gravity  contours  in  Ref.  (18)  are  line-of- 
sight  accelerations  between  the  satellite  and  the  ground  station.  These  line-of- 
sight  accelerations  are  normalized  to  100  km  above  the  mean  lunar  surface. 


Table  I 


Has con  Name 

Lat. 

(Deg. ) 

Long. 

(Deg.) 

Hass 
(xlO-6  Luna 

1. 

Imbrium 

+38 

-18 

20 

2. 

Serenitatis 

+28 

+18 

20 

3. 

Crisium 

+16 

+58 

10 

4. 

.  Nectaris 

-16 

+34 

9 

The  DMC  method  estimates  the  unmodeled  acceleration  vector  which  is  plotted  in 
the  local  spacecraft  coordinate  system  shown  in  Fig.  2.  In  this  figure,  X  is 
the  selenographic  longitude  and  <t>  is  the  latitude.  The  unit  vector  lies 

along  the  selenocentric  position  vector  of  the  satellite,  unit  vector  is 
directed  due  east,  and  is  along  the  spacecraft  latitude  meridian.  The  lo¬ 

cations  of  the  tracking  stations  coincide  with  those  used  in  the  Apollo  missions. 
Generally,  between  four  and  seven  stations  can  observe  the  satellite  simultaneously 
In  all  cases,  the  initial  covariance  matrix,  Pq  ,  was  assumed  to  be  diagonal 

with  the  following  initial  values:  P  =P  =  P  =  o2  =  9  x  lO1*  (meters)2  , 

xx  yy  zz  r 

P. .  =  P. .  =  P. .  =  (a  )2  =  .25  (meters/sec)2,  P  =  P  =  P  =  4  x  10-8 

xx  yy  zz  v  exex  eyey  . 

(meters /sec2) 2  and  P„  ,  =  Pg  g  =  P„  g  =  2.25  x  ICT4  sec2.  The  initial  value  of 

e  is  assumed  to  be  zero.  Finally,  an  observation  interval  of  six  seconds  was  used 
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Numerical  Results 

The  results  obtained  in  the  simulated  study  compare  the  behavior  of  the 
extended  Kalman-Bucy  filter  with  no  model  compensation  with  the  behavior  of  the  DMC 
method  for  both  an  equatorial  orbit  and  an  inclined  orbit.  The  effect  of  the  ob¬ 
servation  accuracy  and  the  unmodeled  acceleration  magnitude  are  considered.  As 
described  in  the  previous  section,  the  true  model  error  is  due  to  the  four  lunar 
surface  mascons . 

1.  Behavior  With  No  Model  Compensation 

For  this  c.'se,  the  observations  were  generated  by  including  the  mascons; 

however,  in  the  estimation  algorithm,  they  were  neglected  in  the  model  used  for 

a  .  The  effects  of  the  unmodeled  accelerations  were  not  accounted  for  through 
m 

the  use  of  the  state  noise  covariance  matrix  Q  in  Eq.  (15).  This  case  then 
represents  the  application  of  the  extended  form  of  the  Kalman-Bucy  filter  with  no 
model  compensation.  The  observation  accuracy  is  1.5  mm/sec.  For  the  retrograde 
equatorial  orbit  (I  =  180°),  the  position  and  velocity  error  norms  are  shown  in 
Fig.  3.  Since  this  is  a  simulation,  the  actual  error  in  the  estimate  is  known 
from 

Ar  =  /(x  -  xT)2. +  (y  -  yT>2  +  (z  -  z^)2 

where  the  (~)  indicates  estimated  values,  and  (  indicates  true  values  used  in 
generating  the  observations.  Fig.  3  also  shows  a  plot  of  the  square  root  of  the 
trace  of  the  first  three  diagonal  el  ;ments  of  the  covariance  matrix  P  ,  i.e.. 


P  +  P 
xx  yy  zz 


Note  that  the  actual  erorr  is  not  bounded  by  the  covariance  matrix.  With  no  model 
compensation,  the  actual  position  error  norm  reaches  a  maximum  of  over  600  meters 
and  the  velocity  error  norm  has  a  peak  at  .5  meters/sec.  The  orbit  with 


I  =  150°  passes  directly  over  two  of  the  simulated  mascons.  With  no  model  com¬ 
pensation,  the  error  norms  are  shown  in  Fig.  4.  The  position  error  reaches  a  value 
of  10,000  meters.  The  estimation  procedure  has  diverged  for  this  case. 

The  position  and  velocity  error  norms  are  shown  in  lieu  of  observation 
residuals.  The  residuals  show  the  effect  of  a  systematic  error  and  will  have  a 
pseudo-periodic  behavior.  Due  to  the  number  of  stations  involved,  the  residuals 
would  require  several  plots  whereas  the  error  norm  requires  only  two  plots. 

Finally,  it  should  be  noted  that  application  of  the  batch  processor  to  these  two 
orbits  will  yield  residuals  of  comparable  magnitude  and  pseudo-periodic  character. 
See,  for  example,  the  results  given  in  Ref.  (14). 

2.  Behavior  Using  the  DMC' Method 

Application  of  the  DMC  method  to  the  I  =  180°  case  yields  the  results  shown  in 
Fig.  5.  Figure  5a  shows  the  radial- component  of  the  unmodeled  acceleration  com¬ 
pared  with  the  estimate  provided  by  the  DMC  algorithm.  The  smooth  line  is  the 
true  acceleration  and  the  uneven  line  is  the  estimate  of  the  true  acceleration. 
Figures  5b  and  5c  show  the  and  components.  The  acceleration  unit  used  in 

*  i 

these  plots  is  the  milligal  (mgal).  Note  that  the  DMC  method  provides  a  very  good 
estimate  of  the  unmodeled  acceleration.  The  position  error  norm  is  illustrated 
in  Fig.  5d  and  the  velocity  error  norm  in  Fig.  5e.  It  can  be  seen  from  Fig.  5d 
that  the  smooth  line,  the  square  root  of  the  trace  of  the  first  three  elements  in 
the  covariance  matrix,  P,  represents  a  conservative  estimate  of  the  position  error 
and  that  the  estimate  of  the  state  of  the  satellite  is  improved  considerably 
over  the  estimate  obtained  in  the  preceding  case,  i.e.,  with  no  model  compensation. 
The  observation  noise  standard  deviation  used  in  both  the  simulation  and  estimation 
computations  was  1.5  mra/sec  as  in  the  preceding  case. 
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For  the  I  =  150°  case,  the  acceleration  estimates  are  shown  in  Fig.  6a  to 
6c.  Since  the  satellite  passes  almost  directly  over  the  location  of  two  of  the 
simulated  inascons,  the  unmodeled  acceleration  is  almost  100  times  larger  than  for 
the  equatorial  orbit.  The  acceleration  estimates  again  yield  accurate  approxima¬ 
tions  of  the  true  acceleration.  Furthermore,  it  can  be  seen  from  Fig.  6d  that  the 
covariance  matrix  provides  a  conservative  estimate  of  the  actual  error.  In  Fig.6e, 
the  true  velocity  error  no™  is  scattered  about  the  covariance  matrix,  thus  the 
covariance  is  indicative  of  the  actual  velocity  error.  The  observation  standard 
deviation  for  this  case  was  1.5  mm/sec. 

3.  Effect  of  Improved  Observation  Accuracy  on  the  DMC  Method 
To  evaluate  the  influence  of  observation  accuracy,  the  case,  I  =  150°,  was 
recomputed  with  an  observation  error  standard  deviation  of  0.15  mm/sec,  i.e.,  an 
improvement  over  the  previous  case  of  an  order  of  magnitude.  The  results  are 
given  in  Fig.  7.  In  Figs.  7a  through  7c,  a  dramatic  improvement  in  the  accuracy 
of  the  estimate  of  the  unmodeled  acceleration  can  be  seen.  Furtheraore,  as  one 

would  expect,  an  improvement  in  the  estimate  of  the  position  and  velocity  of  the 

•  . 

satellite  also  occurs.  Similar  behavior  occurs  for  the  I  =  180°  orbit. 


B*WfiWa«K». ~_ _ . _ _ _ _ 
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Summary  and  Conclusions 

The  results  of  the  preceding  section  demonstrate  that  the  DMC  algorithm 
can  provide  more  accurate  estinates  of  the  state  than  can  be  obtained  with  the 
non- compensated  algorithms.  Furthermore,  the  method  provides  estimates  of  the 
unmodeled  acceleration  time-history  which  are  representative  of  the  true  un- 
modeled  acceleration.  The  accuracy  of  the  acceleration  estimate  depends  on  the 
magnitude  of  the  unmodeled  acceleration,  the  accuracy  of  the  observations  and 
density  of  the  observation  data  set,  i.e.,  the  number  of  tracking  stations  ob¬ 
serving  the  satellite. 

For  a  sufficiently  large  value  of  the  ratio  of  the  magnitude  of  the  un¬ 
modeled  accelerations  to  the  accuracy  of  the  observations  and  for  a  sufficiently 
dense  observation  data  set,  a  precise  estimate  of  the  components  of  the  unr.odeled 
accelerations  can  be  obtained  by  the  Dynamic  Model  Compensation  algorithm  described 
in  this  investigation.  The  estimates  of  the  components  cf  the  unmodeled  accelera¬ 
tions  can  be  used  to  improve  the  knowledge  of  the  lunar  gravitational  potential 
representation. 


i3&l.  • 
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FIGURE  /.  SIMULATED  LUNAR  SATELLITE  GROUND  TRACK 
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Figure  5.  Continued 


Component 


Position  Error  Norm  e.  Velocity  Error  Norm 


Figure  7.  Ia  150°  Unmodeled  Acceleration 
a.  -  0.15  mm/sec 


